{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Geography\n",
    "\n",
    "Datashader contains a `geo` module which contains helper functions which should be familiar to the geospatial community. \n",
    "\n",
    "Some of the functions available include\n",
    "* [Generate Terrain](#ds.geo---generate-terrain)\n",
    "* [Hillshade](#ds.geo---hillshade-function)\n",
    "* [Slope](#ds.geo---slope-function)\n",
    "* [Aspect](#ds.geo---aspect-function)\n",
    "* [Bump](#ds.geo---bump-function)\n",
    "* [NDVI](#ds.geo---ndvi-function)\n",
    "* [Mean](#ds.geo---mean-function)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.transfer_functions import shade, stack\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Terrain Data\n",
    "\n",
    "To demonstrate using these functions, let's generate some fake terrain..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader import Canvas\n",
    "from datashader.geo import generate_terrain\n",
    "from datashader.colors import Elevation\n",
    "\n",
    "W = 1000\n",
    "H = 750\n",
    "\n",
    "canvas = Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))\n",
    "terrain = generate_terrain(canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grayscale value above shows the elevation linearly in intensity (with the large black areas indicating low elevation), but it will look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shade(terrain, cmap=Elevation, how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hillshade\n",
    "\n",
    "[Hillshade](https://en.wikipedia.org/wiki/Terrain_cartography) is a technique used to visualize terrain as shaded relief, illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which is based on slope and aspect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.geo import hillshade\n",
    "\n",
    "illuminated = hillshade(terrain)\n",
    "\n",
    "shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can combine hillshading with elevation colormapping to indicate terrain types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack(shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear'),\n",
    "      shade(terrain, cmap=Elevation, how='linear', alpha=128))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slope\n",
    "[Slope](https://en.wikipedia.org/wiki/Slope) is the inclination of a surface. \n",
    "In geography, *slope* is amount of change in elevation of a terrain regarding its surroundings.\n",
    "\n",
    "Datashader's slope function returns slope in degrees.  Below we highlight areas at risk for avalanche by looking at [slopes around 38 degrees](http://wenatcheeoutdoors.org/2016/04/07/avalanche-abcs-for-snowshoers/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.geo import slope\n",
    "\n",
    "avalanche_slope_risk = slope(terrain)\n",
    "avalanche_slope_risk.data = np.where(np.logical_and(avalanche_slope_risk.data > 25, \n",
    "                                     avalanche_slope_risk.data < 50),\n",
    "                                     1, np.nan)\n",
    "\n",
    "stack(\n",
    "    shade(terrain, cmap=['black', 'white'], how='linear'),\n",
    "    shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "    shade(avalanche_slope_risk, cmap='red', alpha=100), \n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aspect\n",
    "\n",
    "[Aspect](https://en.wikipedia.org/wiki/Aspect_(geography)) is the orientation of slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.\n",
    "\n",
    "Below, we look to find slopes which face close to north."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.geo import aspect\n",
    "\n",
    "north_faces = aspect(terrain)\n",
    "north_faces.data = np.where(np.logical_or(north_faces.data > 350 ,\n",
    "                                          north_faces.data < 10), 1, np.nan)\n",
    "stack(\n",
    "    shade(terrain, cmap=['black', 'white'], how='linear'),\n",
    "    shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "    shade(north_faces, cmap=['aqua'], alpha=50), \n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NDVI\n",
    "\n",
    "The Normalized Difference Vegetation Index (NDVI) quantifies vegetation by measuring the difference between near-infrared (which vegetation strongly reflects) and red light (which vegetation absorbs).\n",
    "\n",
    "For example, when you have negative values, it’s highly likely that it’s water. On the other hand, if you have a NDVI value close to +1, there’s a high possibility that it’s dense green leaves.\n",
    "But when NDVI is close to zero, there isn’t green leaves and it could even be an urbanized area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of *NDVI* ranges from [-1,+1], where `-1` means more \"Red\" radiation while `+1` means more \"NIR\" radiation.\n",
    "\n",
    "Below, we simulate the red and near-infrared bands using `datashader.perlin` random noise with different seeds and frequencies.  Green areas should be those > 0, where higher NDVI values would indicate vegetation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.geo import ndvi\n",
    "from datashader.geo import perlin\n",
    "\n",
    "near_infrared_band = perlin(W, H, freq=(4, 3), seed=1)\n",
    "red_band = perlin(W, H, freq=(32, 32), seed=2)\n",
    "vegetation_index = ndvi(near_infrared_band, red_band)\n",
    "shade(vegetation_index, cmap=['purple','black','green'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bump\n",
    "Bump mapping is a cartographic technique often used to create the appearance of trees or other land features.\n",
    "\n",
    "The `datashader.bump` will produce a bump aggregate that can then used to add detail to the terrain.  In this case, I will pretend the bumps are trees and shade them with green."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.geo import bump\n",
    "from functools import partial\n",
    "\n",
    "def tree_heights(locations, min_val, max_val, height):\n",
    "    out = np.zeros(len(locations))\n",
    "    for i, (x, y) in enumerate(locations):\n",
    "        val = terrain.data[y, x]\n",
    "        if val > min_val and val < max_val:\n",
    "            out[i] = height\n",
    "        else:\n",
    "            out[i] = 0\n",
    "    return out\n",
    "\n",
    "TREE_COUNT = 200000\n",
    "\n",
    "trees = bump(W, H, count=TREE_COUNT // 3,\n",
    "             height_func=partial(tree_heights, min_val=50, max_val=500, height=10))\n",
    "\n",
    "trees += bump(W, H, count=TREE_COUNT,\n",
    "             height_func=partial(tree_heights, min_val=500, max_val=2000, height=20))\n",
    "\n",
    "trees += bump(W, H, count=TREE_COUNT // 3,\n",
    "             height_func=partial(tree_heights, min_val=2000, max_val=3000, height=10))\n",
    "\n",
    "tree_colorize = trees.copy()\n",
    "tree_colorize.data[tree_colorize.data == 0] = np.nan\n",
    "\n",
    "stack(shade(terrain + trees, cmap=['black', 'white'], how='linear'),\n",
    "      shade(hillshade(terrain + trees), cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "      shade(tree_colorize, cmap='limegreen', how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean\n",
    "The `datashader.mean` function will smooth a given aggregate by using a 3x3 mean convolution filter. Optional parameters include `passes`, which is used to run the mean filter multiple times, and also `excludes` which are values that will not be modified by the mean filter.\n",
    "\n",
    "Just for fun, let's add a coastal vignette to give out terrain scene a bit more character. Notice the water below now has a nice coastal gradient which adds some realism to our scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.geo import mean\n",
    "\n",
    "LAND_CONSTANT = 50.\n",
    "\n",
    "water = terrain.copy()\n",
    "water.data = np.where(water.data > 0, LAND_CONSTANT, 0)\n",
    "water = mean(water, passes=10, excludes=[LAND_CONSTANT])\n",
    "water.data[water.data == LAND_CONSTANT] = np.nan\n",
    "\n",
    "stack(\n",
    "    shade(terrain, cmap=['black', 'white'], how='linear'),\n",
    "    shade(water, cmap=['aqua','white']),\n",
    "    shade(hillshade(terrain), cmap=['black', 'white'], alpha=128, how='linear'),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "We've now seen a bunch of datashader's `geo` helper functions.\n",
    "\n",
    "Let's make our final archipelago scene by stacking `terrain`, `water`, `hillshade`, and `tree_highlights` together into one output image: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack(shade(terrain + trees, cmap=Elevation, how='linear'),\n",
    "      shade(water, cmap=['aqua','white']),\n",
    "      shade(hillshade(terrain + trees), cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "      shade(tree_colorize, cmap='limegreen', how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "- Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), pp 406\n",
    "- Making Maps with Noise Functions: https://www.redblobgames.com/maps/terrain-from-noise/\n",
    "- How Aspect Works: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-aspect-works.htm#ESRI_SECTION1_4198691F8852475A9F4BC71246579FAA"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
